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I. Introduction 


The NAVSTAR Global Positioning System (GPS) is a satellite navigation 
system currently under development by the Department of Defense. A 
Joint-Service program office was established in 1974 at the Air Force Space 
and Missile Systems Organization in Los Angeles, California. Tests conducted 
between 1974 and the present have verified the system concept and the 
predicted system performance. GPS will consist of twenty four satellites, 
inclined 60 degrees to the equator, in twelve hour orbits. The GPS signal 
is a spread spectrum signal with two pseudo-random noise codes in quadature. 

A 50 bit per second data bit stream is modulated on both codes. The data 
stream contains system time, satellite clock characteristics data, satellite 
ephemerides, plus other status indicator data. 

The National Aeronautics and Space Administration (NASA), Goddard 
Space Flight Center, has developed a concept utilizing a geostationary 
reference satellite (REFSAT) that broadcasts every four seconds updated 
GPS satellite coordinates. This procedure reduces the complexity of the 
GPS receiver. For an overview of the NASA concept, reference (1) is an 
excellent paper. This paper should be reviewed to understand the REFSAT 
and GPS receiver interface. 

NASA initial direction was to: 

1. Quantify the economic and performance payoffs associated with 
replacing maritime shipborne navigation systems with NAVSTAR, and 

2. Evaluate the use of NAVSTAR for measurements of ocean currents in 
the broad ocean areas of the world. 



- 2 - 


Dun'ng this contract period, a single channel GPS receiver was developed 
by Systematics General Corporation. This unit includes the RF/IF hardware, 
L-band synchronizer and frequency reference, code tracking loop, carrier 
tracking loop and digital processor. The software control and signal 
processing software were developed by Howard University, Washington, D.C. 
There was a requirement to develop a dynamic software program to position 
a ship at sea. Colorado State University (CSU) was redirected to develop 
this software and integrate the software into the REFSAT receiver. 

2. Initial Research 

Cargo and transport ship velocities in the range of 10 to 20 knots can 
have their speed through the water reduced by 10-20% by ocean currents. 

The general trend in ocean currents is quite well understood and has been 
the subject of numerous analytical efforts (2) and also many empirical 
studies with resulting sea charts published by the United States Navy. The 
Pilot Charts include recommended sea-lanes and make the distinction between 
full and low-powered ships(3). 

These currents are time-varying. In particular, both the Gulf Stream 
and the Japanese Current, two of the predominant sources of general oceanic 
circulation for the Northern Hemisphere, can have marked variations, 
particularly along the edges of the mainstream of the current itself. As a 
result, it is worthwhile to seek minimum-fuel ocean passages which rely on 
an understanding of the fine grain, time-varying structure of those ocean 
currents . 

Receivers using the GPS signal have the capability to accurately 
measure ocean currents. This information could be relayed via satellite to a 
central data reduction center where the real-time current data could be 
integrated into the computer data banks. The updated data file could then 
be used to project the minimum sailing time route back to the ship. 
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If the use of GPS would allow an improvement of half a knot in forward 
speed, there is a rather remarkable cost savings in dol 1 ars-per-year for the 
tanker fleet. From Table 1, one can see that a tanker in the 100,000-ton 
class can realize savings of as much as $90,000 per year, if the knowledge of 
ocean currents allows a speed improvement of .46 knots. Figure 1 is an 
extract from the pilot chart of the North Atlantic for July 1978. It shows 
a path which is prescribed for low-powered ships to gulf ports; the path passes 
through regions of extremely confused currents. 


TANKER PERFORMANCE CHART (4) 


Tanker 

Displacement 

Avg. Speed 
Knots 

TTU/Year* 

$/TTU 

3% Speed 
Improvements 

ATTU** 
A3% Speed 

Savi ngs/ 
Year 

19,000 

(FOMALHAUT) 

14.5 

873 

$1 ,276 

.435 kts. 

2.86% 

31 ,858 

52,000 

(VIRGO) 

16.25 

2882 

$ 892 

.487 kts. 

2.80% 

71 ,980 

94,090 

(DRACO) 

15.3 

4810 

$ 666 

.459 kts. 

2.81% 

90,017 


* A TTU is 10^ ton miles 
** 10,000 mile round voyage 

Table 1. Cost Savings Due to Speed Improvements 

Note that the dashed lines on Figure 1 reflect current uncertainty. The 
pilot charts reflect many examples of such unpredictable currents. An 
important case is that in the Gulf of Alaska, in which the area has extensive 
dashed currents for the month of July, 1978. This is the area currently being 
traversed by the Alaskan tanker fleet. Shown on Table 1 are the overall 
savings in dollars-per-year versus speed improvements for three classes of 
ship. 
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The initial research was directed toward constructing an algorithm that 
would give the minimum sailing time route between two points on the earth. 

To specify a position on the earth's surface we need two angles, 6 and 4) . 
Let (0 q,{|)q) be the position of the point of origin and (0^,c()i) the destination. 
Using 0 as a parameter, a function 4>(0) is determined which provides the 
minimum sailing time. 

First an integral is developed which provides the sailing time in terms 
of ocean currents, ship speed and #(6). To find the function that minimizes 
this integral we apply calculus of variations methods to a second order differ- 
ential equation - Euler's equation with boundary conditions. The Euler's 
equation is quite long and nonlinear. Solving the equation numerically, the 
interval from 6 q to 0^ is divided into n subintervals. The first and 
second derivatives of $(0) at each end point of the subinterval are approxi- 
mated by using Taylor's theorem. In place of the Euler's equation we now 
have a system of n-1 nonlinear algebraic equations in the unknowns 
O(0Q+iA0)(i=l ,. . . ,n-l ) plus the boundary conditions. 

Let £ be this system of equations and £ n-1 vector of values of 

$(0Q+iA0). Next the system £(_X)=£ is solved. This system can be rewritten 

£(X) + X = G(X) = X. 

This system is solved by iterating on an initial guess using the equation 

where k = 0,1,3... 

CSU was in the process of implementing the system on the computer when 
we were redirected to develop a navigation filtering algorithm to support the 
REFSAT approach to low cost GPS terminals. 

Some analytical approaches to solving this optimization problem are 


described in Appendix A. 
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3. Development of Navigational Filtering Algorithms . 

CSU was directed to support the REFSAT development by developing naviga- 
tion filtering algorithms. This research includes coordination with Howard 
University and Systematics General. The hardware developed by Systematics 
General Corporation was to be integrated with the digital processing and 
operational software developed by Howard University. Development problems 
delayed the integration. CSU was requested to have the software available 
for integration in the summer of 1981. Unfortunately , funding to the hardware 
and software activities was not available. Funds set aside for the integra- 
tion of the navigation software have been returned by CSU to NASA. 

The navigation software was to be resident in the REFSAT receiver 
processor. As the digital processing and operational software developed, it 
appeared that insufficient space would be available, both from a throughput 
and a memory point of view. It was therefore concluded that the software 
should not be developed independent of the operating REFSAT receiver. 

Information on the Navigation Function Driver (NAVFUNC) has been received 
from Howard University. It is possible that some of the data calculated in 
the subroutines could be used to reduce the computation time of the navigation 
filter. If the NAVFUNC data is not retrievable, then a method similar to the 
one used by A. J. Van Dierendonck (5) must be used to obtain the NAVSTAR 
satellite coordinates. 

Time to calcualte the navigation solution is very limited. It is estimated 
that it will take from 5-10 seconds to calculate the necessary variables 
needed for the navigation algorithm. The solution of the ship's position will 
take from 15-30 seconds. For example, to solve for the linear system repre- 
sented by 
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Hx = y (rank 5) 

takes from 1,7 to 2.1 seconds to solve. A listing of the program to solve this 
system is listed in Appendix B. The dead reckoning portion will take another 
5-10 seconds, for a total computational time of slightly less than one minute. 
Times are based onan 8080/Z80 based microprocessor 2MHZ clock). 

The accuracy of the position and velocity of the ship will be impacted 
by the ship's motion: headway, surge, leeway, sway, heave, roll, pitch and 

yaw. Of these motions, roll will have the greatest adverse impact. Large 
ships (over 150 meters in length) can have roll periods up to 13 seconds. 

By taking a sampling period much greater ('^> 5x) than the period of roll the 
short term vessel's motion could be treated as increased noise. 

The initial navigation algorithm will be based on the work of Noe, Myers 
and Wu (6), The navigation solution will be strongly dependent upon system 
noise. Since the REFSAT receiver performance has not been fully defined, 
the navigation accuracy can not yet be determined. 

4 , Future Efforts 

The incorporation of the navigation filter into the REFSAT receiver should 
take approximately one month. If the receiver is not operational or persons 
familiar with the receiver are not available for consultation then the period 
of effort would extend several months. CSU has the personnel capable of making 
the receiver operational and would be pleased to continue this effort if and 
when funding becomes available. 
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Appendix A 


The trajectory optimization problem discussed in Section 2 is described 
in more detail in this appendix. A simplified analytical formulation is 
presented based on a kinematical analysis and a closed form solution to this 
simplified problem is presented. Next, various complications are added and 
the form of the optimal controller for each variation is presented. Finally, 
a dynamical formulation is presented and its solution is discussed. 

The problem of determining the optimal trajectory for a ship may be 
described in simplified terms as follows. Consider the situation represented 
on Figure A.l. The problem is to determine the steering angle, 6(t), that 
results in a transfer of the ship from its initial position (L^,a^) to 
a final position (L^,a^) in minimum time. Without loss of generality we may 
assume that the final latitude, L^, is zero. Also, without loss of generality, 
we may assume that the motion takes place on a unit sphere and that the 
ship's speed relative to the water is unity. 

Model No. 1. No Ocean Currents 

For this initial problem we assume that there are no ocean currents. 

This provides a starting point for which we can obtain a closed form analytical 
solution. We may then add the ocean currents to the model and attempt to 
build a solution in steps. 

The kinematic equations may be written as 



z 


ORlGiNAL 23 

OF POOR QUALITY 



Figure A.1 . Ship's Trajectory 
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^ = sin 6 


da _ cos 6 
dt ~ cos L 


The payoff function may be defined as 


4> = -t, 


(A.l) 
(A. 2) 


(A. 3) 


where is the final time. We wish to determine 0(t) such that cj) is 
maximized as this is equivalent to minimizing the time required for the 
transfer to take place. The terminal constraints are 

= L(t^) = 0 (A. 4) 


\l>2 = a(t^) - = 0 


(A. 5) 


The Hamiltonian function is 


H = A, sin 9 + A 


cos 6 
a cos L 


(A. 6) 


A necessary condition for optimality is that the control, 9(t), maximizes the 
Hamiltonian. Therefore the optimal control is determined unambiguously by 
the equations 

9 = \ A? + A^ / cos^T (A. 7) 

La 

cos e = (X_^/cos L) / '' cos' L (A. 8) 
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Equations (A, 7) and (A. 8) specify the optimal control in terms of the adjoint 
variables and These, in turn, are determined from the following 
equations . 


^ 9H 
dt ‘ 9L 

2 

= - A cos 6 sin L / cos L 
a 


dA 
a 

dt 


3H 

9L 


= 0 


(A. 9) 


(A. 10) 


The adjoint variables are determined at the final time, t^, from the trans- 
versal ity condition 



where 


T 

X 


(L,a). From equation (A. 3), 


(A. 11) 


3X 


and from equations (A. 4 ) and (A. 5) 


(0 


0 ) 



■ 1 

o' 

^f 

0 

1 


(A. 12) 


9i{; 

9X 


(A. 13) 
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Therefore 


U|_. = (-V, -v^) (A. 14) 

T "f 

The components of v represent the sensitivities in the payoff function, (j), 
to changes in the components of the state vector at the final time. 

Consequently 


v-| = -sin 6^ 


and 


= -cos 


and therefore 


Aj^(t^) = sin 0^ 


(A. 15) 


and 


(t^) = cos 6^ 


(A. 16) 


An additional necessary condition is that 


dt 


[ <P - 


] 


= 0 


(A. 17) 


Expanding, 


^9x ax^ dt '3t 


'' w ^ 


= 0 


(A. 18) 
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Recalling equations (A.1), (A. 2), (A. 6), and (A. 11), and using the fact that 
= -1 and 4+ = 0> this reduces to 


H|t^ = 1 (A. 19) 

Also, since the Hamiltonian is not an explicit function of time it must be 
constant along an optimal trajectory. Therefore, necessarily 

H E 1 

along the optimal path. This provides us with the following equation for 

= 1 - cos^ 0^ / cos^ L (A. 20) 


The above system of equations may be solved in closed form; the solu- 
tions are 


L(t) = -sin“^ [sin 0^ sin(t^ - t)] 


a(t) = 01 ^ - tan~^ [cos 0^ tan (t^ - t)] 


\(‘> ' ^0= (V - ' W + ootHtf - t) 


A^(t) = cos 0^ 


(A. 21) 
(A. 22) 
(A. 23) 
(A. 24) 


and 


0(t) = tan 


[tan 0^ cos (t^ - t)] 


(A. 25) 
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These solutions may be verified by substituting them back into the differential 
equations. For any given value for a^, the above equations define a two- 
parameter family of trajectories on the sphere. The two parameters are 
0^ and t^. 

Given any arbitrary state (L,a) at some time t < t^ one may solve for 
0^ and t^ in terms of L, a, and t. The results are 

Q - 4 . - 1 r -tan L n 

f (A. 26) 

and 

t^ = t + cos'^ [cos L cos (a^ - a)]. (A. 27) 

These relations for 0^ and t^ may be substituted into equations (A. 23), 

(A. 24) and (A. 25) to obtain the optimal values for the adjoint variables and 
the control. The results are 


Al* (L, 


= - si" L / L + tan^ (a^ - o) 


(ft. 28) 


y (L.a.t) = cos L tan (a^ - a) / '' s,„2 l + tan" (a,, - o) 


(A. 29) 


and 


e* (L.o.t) = tan-' [ 3 

The optimal value for the payoff function is 


-1 


cp* (L,a,t) = -t^* (L,a,t) - cos' [cos L cos (a.p - a)] 


(A. 30) 


(A. 31 ) 
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Model No. 2. Constant Current (East or West) 

We next consider a situation in which there is a constant ocean current 
in the easterly or westerly direction. In this case the kinematical equations 
are 


dt 


sin 6 


(A. 32) 


da _ cos 6 ± k 
dt cos L 


(A. 33) 


The Hamiltonian is 


H = X, cos 9 + X 
L a 


f Cos 6 ± k x 
^ cos L ' 


(A. 34) 


The Hamiltonian is again maximized by the control 

X, X cos L 

sin e = - — ■ , cos 0 = ■ ■ ^ (A. 35) 

*^X| ^ + X ^ / cos^ L X, ^ + X ^ / cos^ L 

La La 

That is, the optimal control is determined by the same equations as for the 
first problem. 

Also, we still have that 


dt 


T^ 


dt t 


= 0 


(A. 36) 




- V 


Ml ^ 

3X^ dt |t 


ht ■ 


dip 

9t 




= 0 


(A. 37) 
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Therefore 




_ T 


= d 


(4.38) 


But 


^4> - 1 M = n 
3t ’ 3t 


(4.39) 


Also 


F k = H 


Therefore 




(A. 40) 


• dH 

Also H ^ H(t) ; •• 0 along the optimal trajectory; also, we have that 

H E 1 along the optimal trajectory. Hence 


X| sine + X ~-T + ^ ^ 

L a cos L — a cos L 


(A. 41) 


The adjoint equations are 


dX, 


L -3H 


dt 9L 


- X^(cos e + k) 

cos L 


(A. 42) 


and 




(A. 43) 


Also, we still have 


X (t,) = 
a ' f 


at. 

= cos 0. 

aa. f 


(A. 44) 
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A1 so, from H = 1 , 


^ + A ^ / cos^ L 
L a 


(X^ / cos L)^ 

X, ^ + X ^ / cos^r 

L a 


X k 
- cos L 


1 



X 


L 


2 


+ 


Xg^ 

cos^ L 


, , 2 X k 

+ X , g 

c^ - L 


= 1 


_ X k o 

(1 + — = 1 + 
cos L 


X^k^ 
g 

cos^L 


a + 


2X k 
g 

cos L 


{A. 45) 


The solution proceeds along the same lines as for the first problem. 
Model Number 3. Constant Current (North and South) 

Let = constant in the north or south direction. Then 


^ = sin 9 + k 

dg _ Cos 9 
dt ~ cos L 


H = X, (sin 9 + k) + X 

I — ' g cos L 


= X, sin 9 + r cos 9 + kX, 
L cos L — L 


(A. 46) 
(A. 47) 


(A. 48) 


Again the control law will remain the same. 
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That IS 


sin 6 


2 + Xh 


, cos 0 


cos^ L 


\ / cos L 


(A. 49) 


cos^ L 


But the adjoint equations are 


^ _ 3H 
dt ~ ■ 3L 


A cos 0 
■a . , 

2 L 

cos*^ L 


(A. 50) 


dA 
a 

dt 



(A. 51 ) 


Again we see that the addition of a constant current to the model results in 
the same form of solution. Clearly the second and third models may be 
combined. 

Model No. 4. Linearly Increasing Currents (East and West) 

Now suppose that the ocean current is a linearly increasing function of 
latitude. That is 


= + kL (A. 52) 

Then the system equations are 

^ = sin 0 (A. 53) 


da _ cos 9 t kL 
dt ” cos L 


(A. 54) 
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The Hamiltonian is 


H = X, sin 6 + X (cos 6 + kL)/cos L (A. 55) 

L a ~- 

The optimal control law still remains the same as before; the adjoint variables, 
which actually determine the optimal control as a function of time are given 
in this case by 


dX, I kX 

= - A (cos 0 t kL) i 

dt a' - ' ^^^2 cos 


COS L 


dX 
( 

dt 


a 


= 0 


The solution technique proceeds in the same manner as before. 

Model No. 5. Linearly Increasing Currents (North or South) 

In this case consider 


(A. 56) 


(A. 57) 


± ka 


(A. 58) 


Then the system equations are 


dt 


:rr = sin 0 + ka 


(A. 59) 


da _ cos 6 
dt cos L 


(A. 60) 


The Hamiltonian is 


H = X. (sin 0 + ka) + X cos 0/cos L 
L — a 


(A. 61) 
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Again there is no change in the optimal control law.' The adjoint equations are 


and 



p 

A cos 6 sin L/cos i 
a 



+ kA|^ 


(A. 62) 


(A. 63) 


Thus we are no longer have that A^ is constant along an optimal trajectory. 


Model No. 6. Dynamical Model 

The problem of determining the optimal steering angle to minimize transfer 
time is similar to the problem of steering a ship to minimize propulsion losses. 
This latter problem has been addressed [A.l], [A. 2] and is briefly discussed 
below. Comsider the representation depicted in Figure A. 2. The ship's 
equations of motion may be written as 


m - 

r) - r. 

V r 

- r 

p 

0 

N‘ 

V 

(I - N') 
zz r 

- (N' t 
P 

I )0 
xz 

K' 

V 

*(K’ + I ) 
r xz 

(I 

XX 

K' )0 
P 

0 

0 

0 

i 


where 


r 

P 


'y Y - mU Y Y^' 


V 




Y 

V r p $ 




6r 


D 

N N N N. 


r 




N 

V r P t 




6r 


D 

Kr ^ K K. 

V r p (P 


P 

+ 


6r 

"d 

j 

o 

O 

O 




_0 


_0 


(A. 64) 


\p = yaw 
V = sway 
p = roll 

and Yq, Np, Kp represent the external disturbances on the ship due to the 
seaway. In addition, the ship's lateral motion is defined by 


y = V + Uij; 


= r. 


(A. 65) 
(A. 66) 
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Fig. A. 2 Ship coordinate system-body axes coordinates, 
forces, and moments. Abstracted from [A.l]. 

Equations (A. 64), (A. 65), and (A. 66) may be written in state space representa- 
tion as 

AX = BX + Cu + DW , (A. 67) 

or, rewriting in standard form 

X = FX + B^u + GW (A. 68) 

where F=A'^B, B^=A"^C, and G=A"^D. The state vector is X, the control vector 
is u, and the disturbance vector is W. 

Representative open loop lateral plane system eigenvalues for a highspeed 
containership at full-load design condition for speeds of 16, 23, and 32 
knots are [A.l]: 



16 Knots 


0.00 


ORIGINAL PAGE IS 
OF POOR QUALITY 


-0.0584 X 10"^ 
-0.792 X 10"^ 

23 Knots 0.00 

-0.839 X 10'^ 
-0.114 

32 Knots 0.00 

-0.117 X 10"^ 
-0.158 


A simulation model for this system has been developed [A.l] and is represented 
pictorial ly by the block diagram shown on Figure A. 3. 



Figure A. 3 Nonlinear simulation model. Abstracted from [A.l]. 


The problem is to determine the controller that will minimize propulsion 
losses due to steering. These losses may be shown to result from excess 
power consumption per unit distance travelled in the X-direction caused by 
the added resistance due to steering. This is related to the instantaneous 
surge due to steering, which is 
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AX = (m + (p/2) LAX;^) vr + (1/2)](p/2) AX;^] / (A. 69) 

+ (l/2)[(p/2) AX'g^ U^] 6^ 

where m denetes mass pf ship; p, density of sea water; L, ship length between 

perpendiculars; A=L ; U, ship's water speed; v, sway velocity of ship; r, yaw 

rate of ship; 6, rudder angle; X^^, force coefficient due to yaw/sway (+ve); 

X'^o, force coefficient due to rudder angle (-ve); X' , force coefficient due 
00 vv 

to sway. From this, the mean surge relevant to steering may be written as; 


AX = [m + (p/2) LAX^^] (v^ r^/2) cos ((J)^ - (J)^) 

+ (1/2) [(p/2) AX' ](v?/2) 

V V a 

+ (1/2) [(p/2) AX^^ U^](6g/2) (A. 70) 


Since the ship motions resulting from seaway disturbances are oscil latory, a 
performance criterion for added resistance due to steering may be formulated 
as: 


J = (1/2) / l-X" vr + nv^ + 6^] dt 

0 


where 


X" = 2[m + (p/2) LAX;^]/[(p/2) AX'^ U^] 


(A. 71) 


(A. 72) 
(A. 73) 


both inversely proportional to speed squared. 

The simulation program described in Ref. [A.l] has been used for the 
purpose of evaluating controllers to minimize the cost function, J, defined 
by equation (A.71 ) . 
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By combining the two methods of analysis described in this appendix the 
global optimal controller may be determined. 
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APPENDIX B 


This appendix includes the listing of the program to solve the linear 
system Hx = y of rank five. Two routines are used to solve the same linear 
system. 

Subroutine SOLVED is written in the normal FORTRAN manner using Do-loops, 
double subscripting, etc. Subroutine SOLVE is written in a more primitive 
way. SOLVE evaluates all subscripts by single additions, uses no Do-loops, and 
isolates floating point operations in a way that will facilitate an easy 
transition to an assembly language version of the code. 

Code developed for the navigation filter will incorporate the features of 
SOLVE. If the time required to obtain a ship's position needs to be reduced 
then the program could be written with a minimum of effort. 

On the MDS-230 system MHz) the approximate times required for the 
two methods are: 

SOLVED 2.1 sec. 


SOLVE 


1 .7 sec. 



C: 

C 
C 
C 
C 

PROGRRM NS 

[:• I M E r J S I 0 N H 2 4 ^ H T H 0 <. 5 .. 6 > .. H T H C 5 .. 6 > .. V C 2 4 .■ X 0 
DRTR IWX ex 
c 

OPEN Uf J I T=^ 1 W.. F I LE = - : FI : UUTPT - > 

RERDkS.. =+^. END=S9y.. ERR=999> ITER 
URITE-:: IN,. ITER 

5 FuRMRTc:x NUMBER OF ITERRTIONS ON ERCH LOOP = 
C 

DO 20 1=1.. 24 

V I ;:■ = SI N F LORT c; I > X 24. ^ 3 . 1416 > > 

DO 10 .J=l.. 5 

H I .. .T = FLORT I +.J > + 1. XFLORT C I +.J > 

10 CONTINUE 
20 CONTINUE 

C 

DO 40 1=1.. 24 

nf:: I TE c; i i-i.. 3e > < h i .. J > .. J =i . 5 > .. v c i ;> 

30 FORMRT«;;5>:;iX.. E10. 4>.. 4X.. E10. 4> 

40 CONTINUE 
C 

PRUSE ■■■ BEGIN FIRST LOOP 

DO 50 1=1.. ITER 

CRLL SOLVE0'-.;H.. HTH0.. V.. XG> 

50 CONT I NUE 

PRUSE ' END OF FIRST LOOP 


>JRITE':; IW.. 60 > 


60 

FORr-IRT 
DO 70 1= 

1.. 5 





WRITE-:: IN.. 

30 

-::HTH0-: 

: I.. j;c .j=i. 6> 


70 

CONT I NUE 
WR I TE I W.. 

80 

X0 i 

-.. 1=1.. 5:;- 


00 

FORMRT X -• 

SOL LIT I ON 

VECTOR X0- X5-:;iX 

.. Ell. 5 


PRUSE BEGIN SECOND LOOP 

DO 90 1=1.. ITER 

CRLL SOLVE tiH.. HTH.. V.. X> 

90 CONTINUE 

PRUSE END OF SECOND LOOP 

C 

DO 100 !t=l.. 5 

WR I TE I W/ 30 > c! HTH C I .. J > .. J=l.. 6 > 

10© CONTINUE 

WR I TE <; I w. 110 > <; X <; I ;c i =i.. 5 > 

110 FORMRTCX' SOLUTION VECTOR X" X5 IX. Ell. 5 > X 
C 

999 STOP 

END 

C 

C 

C 

C 

C 
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SUBROUTINE SOLVEO H.. G. V.- 
D I MENS I Ot J H 24.. 5 > .. G C 5^ 6 > V C 24 > .■ X < 5 > 

DO ±00 J=±.. 5 

DO Sv3 1=1.. 5 
I F I . GT. -J GO T O U0 
S= 0. 

DO 50 K = l.. 24 

S= S + H<K.. .J > =+'H K.‘ I 
50 CONT I NUE 
G <; I . J > = S 
I>= S 

S0 CONTINUE 
100 CONTINUE 

DO 150 1=1.. ?:• 

DO 120 .J = l.. 24 
S= S + J.. I > 

120 CONTINUE 
G f! I .. 6 > = S 
150 CONT I NUE 


DO j;;0O .J = l.. 4 
R0= -1. .■■■Gc;.j.. .j;> 

JP1= .j+i 

DO 280 I = ..TP1.. 5 
P! = P'. 0 '•+■• G I .. ■ J .■* 

DO 250 K=JP1.. 6 

G I .. K = R :+• G '•! J .. + G I .* 

CONT I NUE 
CONT I NUE 
CONT I NUE 

I='i 5 y = G 5.. 8 .■•■'G 5.. 5 
DO 350 11=1.. 4 

1= 5- I I 
IP1= I+l 

S= G I .. 6 > 

DO 320 .J=IP1. 5 

S= S ~ G I .. .J .> '.. -J 

320 CONTINUE 

X'.;i .:'= s. ‘G<;i.. I > 

350 CONTINUE ’ 

RETURN 

END 


250 

280 

300 
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i.ucir- I i.i 4c. U - , 

DIMENSI0r4 Hc;i>.. VCii), X<1> 

INTEGER GRDR. HRDR.- Vr6r.. GPTR. GTBRSE. HBRSE.- GTPTR. 
+ HTBRSE.. HPTR.. HTPTR. VBRSE. VPTR 


integer XRDR.. R0PTR,. GJBRSE, GIBRSE.. RPTR.. 
+ GIPTR,. GJPTR.. GBRSE.. X.PTR 

INTEGER=+1 I.J. K 


HRDR= 1 

grdf-!= !i 

VRDR= 1 
XRL'R= 1 


UF'TF!= GHC'R ~ 1 
GTEiR 1'.-E:= GRD'R ~ t' 
HBRSE= HRDR - 25 

J= € 

10 J= J-1, 

IF<J. EG!. Hi > GO TO iOGi 
GTBRSE= GTBRSE + 1 
GTPTR= GTBRSE 
HTBRSE= HRDR - 25 
HBRSE= HBRSE + 24 

1= 6 

20 1= I-l 

EG!. e> GO TO 10 
GPTR= GPTR + 1 
GTPTR= GTPTR + 5 
HTBRSE= HTBRSE + 24 
IF Cl. LT. J > GO TO 20 


HPTR= HBRSE 
HTPTR= HTBRSE 
S= 0. 

K= 25 
20 K= K-1 

IF<:K. EG!. 01 > GO TO 5Gi 
HPTR= HPTR + 1 
HTPTR= HTPTR + 1 
S= S + H HPTR >:+ H';: HTPTR 
GO TO 30 . 

50 G < GPTR .> = ."S 

I Ft; I. EG!. J;> GO TO 20 
G<: GTPTR :>= S 
GO T O 201 
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c 


100 HPTR= HRDF.: - 1 ORIGINAL PAGE IS 

VE:RSE= VRDR - 1 ’ OF POOR QUALITY 

C 

J= 6 

110 J= J -1 

IFCJ. EQ. 0 > GO TO 140 
GF'TR= GPTR + 1 
VPTF:= VBRSE 
S= 0. 

1= 

12:0 1 = 1-1 

I F '■ 1 . E G!. Gi ' GO TO 13 0 
HPTR= HPTR + 1 
VPTR= VPTR + 1 


S= S + HGHPTR:;':+ V': VPTR> 

GO TGj 12:ti 
1 0 G •: G F‘ T F: = In- 
GGi TO llG-1 
C: 

C: 

14G.i F!0F'TF;= GRC^Fl — €> 

G-JE!RSE= GRDF: — E 
C 

J = 0 

210 J= J-t-1 

1F<;.J EG:i. 5;> GO TO 3:0Gii 
R0PTF:= R0PTR + 6 
R0= -1. /g-;:r0PTR> 

RPTR= R0PTF: 

C 

G-TEiRSE= G.JEiRSE E 
GIE:RSE= G.JBRSE 
C 

1= J 

22£i 1= I-*-l 

IF<;i. EQ. 6> GO TO 21Gi 
RPTR= RPTR + 1 
F: = F: 0 =+= G . F: F' T R y 
GIBRSE= GIBRBE + 1 
GIPTR= GIBRSE 
GJPTR= GJBRSE 
C 

K= J 

23:0 K== F>1 

IFCK. EQ. 7> GO TO 220 
GIPTR= GIPTR + 5 
GJPTF:= GJPTR + 5 

G G G I F’TF: > = F‘>t-G G G JF'TF: GCGIF'TR.'' 

GO TO 23:0 
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30u GEiRi'E= GRDF-! + 30 
1=6 

210 1= I-l 

IFc;i. EQ. 0> GO TO 340 
GBRSE= GBRSE - 1 
GPTR= GBRSE 
XPTR= XRDR + 5 

&= gcgptr;- 

. J = 6 

220 J= J-1 

IF<! J. EG-!. I> GO TO 220 
GPTR= GPTR - 5 
XPTR= XPTR - 1 

s= s - g<;gptr:;-:+^x<kptr> 

GO 1 0 3;'20 

220 >=:ptr= XPTR - 1 
GPTR= GPTR - 
X <; XPTR = S/G C GPTR > 

GO T O 210 

240 RETURN 


EN[;- 


SOLUTION VECTOR 

12262E+00 . 2S»912E+00 62771E-01 2115S'E+00 

■SOLUTION VECTOR X 

12362E+00 . 2-S'5'12E+00 62771E-01 2115‘3E+00 


. 11135E+00 


. 111-35E+0O 



